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ABSTRACT 

This paper presents recent work on developing methods 
for analyzing radiation heat transfer between diffuse- 
gray surfaces using ^-version finite elements. The work 
was motivated by a thermal analysis of a High Speed 
Civil Transport (HSCT) wing structure which showed 
the importance of radiation heat transfer throughout the 
structure. The analysis also showed that refining the 
finite element mesh to accurately capture the 
temperature distribution on the internal structure led to 
very large meshes with unacceptably long execution 
times. 

Traditional methods for calculating surface-to-surface 
radiation are based on assumptions that are not 
appropriate for p-version finite elements. Two methods 
for determining internal radiation heat transfer are 
developed for one and two-dimensional /^-version finite 
elements. In the first method, higher-order elements are 
divided into a number of sub-elements. Traditional 
methods are used to determine radiation heat flux along 
each sub-element and then mapped back to the parent 
element. In the second method, the radiation heat 
transfer equations are numerically integrated over the 
higher-order element. Comparisons with analytical 
solutions show that the integration scheme is generally 
more accurate than the sub-element method. 
Comparison to results from traditional finite elements 
shows that significant reduction in the number of 
elements in the mesh is possible using higher-order (in- 
version) finite elements. 


INTRODUCTION 

Background 

One of NASA’s main goals is to provide the United 
States’ Aeronautics Industry with the technology it 
needs to lead the international aerospace industry into 
the next century. This includes developing enabling 
technologies, one of which is to “provide next- 
generation design tools and experimental aircraft to 
increase design confidence, and cut the development 
cycle time for aircraft in half.”[l] This paper presents 
initial work on one of these tools, methods for accurate 
thermal analysis of aircraft structures. The goal of this 
work is to enable thermal analysis of large-scale 
components or full aerospace vehicles on a mid-level 
workstation. 

NASA sponsored several projects aimed at reducing 
barriers to commercially viable high-speed civil and 
space transportation. The High Speed Research (HSR) 
Program had a goal “of reducing the travel time to the 
Far East and Europe by 50 percent within 20 years, and 
to do so at today’s subsonic ticket prices.”[l] These 
aircraft will experience higher in-flight temperatures 
due to the increased rate of aerodynamic heating 
associated with high-speed flight. The Reusable 
Launch Vehicle (RLV) program’s goal of “reducing the 
payload cost to low-Earth orbit by an additional order 
of magnitude” [1] challenges designers to come up with 
both materials and structural concepts that can 
withstand the reentry environment while minimizing 
the weight of the vehicle. RLV designs are single- 
stage-to-orbit vehicles that use cryogenic propellant. 
The propellant tanks are an integral part of the vehicle 
structure, which, in addition to holding the cryogenic 
fuel, must withstand the elevated temperatures 
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associated with the reentry. In addition, these materials 
must also withstand the thermal and structural cycling 
of repeated flights. Many of the materials that can 
withstand the elevated temperatures of high-speed flight 
do not have high thermal conductivity to efficiently 
dissipate the heat throughout a structure. Thus other 
modes of heat transfer, radiation in particular, become 
more significant. [2] The importance of radiation heat 
transfer as well as the difficulty involved in its 
modeling was demonstrated by Ko et ah [2] in their 
work on the Space Shuttle thermal analysis. These 
issues are directly applicable to the design of high- 
speed commercial aircraft as well as the RLV. 

NASA’s Space Program also has an interest in 
improved methodologies for radiation heat transfer. 
Besides the access to space issue mentioned above, 
radiation is a dominate mode of heat transfer in most 
orbiting space structures and satellites. Chin et al. [3] 
discuss some of the problems associated with thermal 
modeling of spacecraft, and in particular the problems 
associated with computation of radiation heat transfer. . 
One of the problems highlighted in reference 3 is the 
assumption of isothermal, constant radiation heat flux 
surfaces used in most radiation computations. This 
assumption tends to under-predict the temperature 
gradients in the structure. In addition, the assumption 
of constant radiation heat flux over an element may 
degrade the accuracy of the calculated temperature field 
especially where partial shading (blockage) occurs. 

The desire for better tools led to an initial study of the 
heat transfer in a High Speed Civil Transport (HSCT) 
wing. The study confirmed the significance of radiation 
in the heat transfer throughout the vehicle’s wing 
structure. The study also demonstrated the relative 
difficulty of such an analysis, especially when the 
ultimate goal is to compute the temperatures throughout 
the whole wing without reverting to reduced models. A 
summary of this study is given in the following section. 

Case Study: High Speed Civil Transport Wing 
Thermal Analysis 

As part of the HSR Program at NASA Langley 
Research Center (LaRC), a thermal analysis of a HSCT 
wing was undertaken. The purpose of the analysis was 
to determine the capabilities of the methods currently 
used for thermal analysis of aerospace structures and 
determine what, if any, areas of improvement were 
required. 

The wing geometry model used for the thermal analysis 
is shown in Figure 1. The wing is approximately 113 
feet long at the root and 55 feet wide at the trailing 
edge. The wing skin was assumed to be constructed of 


hat-stiffened corrugated panels made from titanium. To 
model the three-dimensional hat-stiffened skin with 
two-dimensional elements, equivalent properties 
(density, capacitance, and thermal conductivity) were 
derived. For simplicity, this construction was also used 
for the internal ribs and spars. 



Figure 1: HSCT wing geometry model. 


A five-hour flight trajectory, representative of a 
commercial airline or transport route, was used in the 
analysis. The majority of the trajectory, 3.8 hours, 
consisted of a mach 2.4 cruise at an altitude between 
60,000 and 70,000 feet. Aerodynamic heating rates 
were determined using LANMIN, NASA Langley’s 
version of the MINIVER computer code which uses 
engineering relations to calculate the aerothermal 
heating (or cooling) to surfaces. The heating rates 
applied to the upper and lower wing surfaces were 
generated by MINIVER for each element in the mesh 
using flat plate boundary layer relations based on the 
element running length and local flow angle. 
MINIVER used inviscid flow relations to obtain the 
undisturbed flow conditions behind the bow shock off 
the fuselage nose and wing leading edge. 

Structural temperatures of the wing throughout the 
trajectory were computed using MacNeal-Schwendler 
Corporation’s P/Thermal 1 [4] software on a Silicon 
Graphics Workstation with an R4000 CPU, 96 Mbytes 
of memory and one Gbyte of disk space. The initial 
mesh consisted of roughly one element for each 
geometry surface shown in Figure 1 resulting in 154 
nodes and 237 elements. Solutions obtained for this 


1 The use of trademarks or names of manufacturers in 
this report is for accurate reporting and does not 
constitute an official endorsement, either expressed or 
implied, of such products or manufacturers by the 
National Aeronautics and Space Administration. 
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mesh both with and without internal radiation showed 
very little difference in temperature distribution. Note 
however that all the nodes in this model are located on 
the wing surfaces (at the intersections of the lines in 
Figure 1) where the heat transfer is dominated by the 
aerodynamic heating. Thus all the computed 
temperatures respond very quickly to changes in the 
aerodynamic heating rates. The lack of nodes along the 
internal structure, where temperatures are not directly 
affected by the aerodynamic heating and thus respond 
slower, results in a misleading temperature distribution 
for the wing structure. 



Figure 2: HSCT wing, refined mesh and 
temperature contours (Temperatures in °F). 

A second analysis was performed with a mesh refined 
to have elements with edge lengths of 18 inches on both 
the internal structure and wing surfaces. This mesh size 
was driven by the desire to have at least two elements 
across the height of the internal structure while keeping 
the overall size of the model within reason. The mesh, 
which consisted of 3624 nodes and 4037 elements, is 
shown in Figure 2 overlaid on a contour plot of the 
lower surface temperature at the beginning of cruise. 
The solution shown in Figure 2 did not include internal 
radiation heat transfer because the view factor 
calculation never ran to completion. It is also 
interesting to note that the view factor computation 
required over 600 Mbytes of disk space before 
crashing! Neglecting internal radiation resulted in a 
peak internal temperature of less than 150°F and a 
temperature gradient of 220°F across the internal 
structure at the beginning of cruise. The large 
temperature difference between the external and 
internal structures raises the possibility of significant 
radiation exchange. 

The problems with the view factor calculation as well 
as the desire to get a more detailed temperature 
distribution for the internal structure led to the 
development of a wing box model. This model allowed 


for a more detailed analysis of a smaller section of the 
wing, in this case a box defined by adjacent ribs and 
spars and the corresponding upper and lower surface 
sections. Solid elements with equivalent properties 
were used for all sides of the box so that temperature 
gradients through the thickness of the panels could be 
determined in addition to the effects of internal 
radiation. The finite element mesh used for the wing 
box thermal analysis is shown in Figure 3. 



Figure 3: Finite element mesh used for the thermal 
analysis of a wing box (Note node 957 is located on 
the lower wing surface). 

The model was first run without internal radiation and 
the transient temperature response of the nodes labeled 
in Figure 3 are shown in Figure 4. Internal radiation 
was added for the next run and the results are shown in 
Figure 5. These figures clearly show the significant 
impact of internal radiation heat transfer on the 
temperatures of the internal wing structure. Ko et ah 
have also shown the importance of internal radiation in 
their work on the Space Shuttle[2]. Their results for the 
in-plane temperature distributions for the upper and 
lower surfaces suggest that higher order basis functions 
might be well suited to this application. This behavior 
coupled with the difficulties encountered using 
traditional methods lead to the idea of applying 
hierarchical /?-version elements to radiation heat 
transfer in enclosures. 
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Figure 4: Transient temperatures for wing box 
model with no internal radiation. 



Figure 5: Transient temperatures for wing box 
model with internal radiation. 

FINITE ELEMENT FORMULATION 

The finite element method and its application to thermal 
problems is well established. For the interested reader, 
Huebner etal. [5] provides a detailed formulation of the 
problem. Here it will simply be noted that the approach 
begins by subdividing the problem region into elements 
and approximating the temperature in each element 
with some function in the form of 

T {e \x,y,z,t)= I N (x,y,z)T (t) (1) 

i = l * 

where are interpolation functions and T l nodal 
quantities related to the temperature at time t. The 
corresponding finite element equations can be written 
(see [5]): 

( 2 ) 
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with r, c p , and k being the material density, specific 
heat, and thermal conductivity, respectively. Boundary 
conditions included in the above definitions are: 
specified heat flux q s on surface S h convection on 
surface S \ with a convection coefficient h and a fluid 
temperature 7}, and net radiation heat flux q r on surface 
S 3 . Q represents any heat source or sink in the material. 

In general, any given element in the solution domain 
will not have all of the above terms. For example, only 
elements on the surface of the domain will have the 
terms associated with the boundary conditions, and 
even then they will only have the terms associated with 
the boundary conditions that apply to the element. The 
terms are calculated individually for each element in the 
domain and then assembled into a global system of 
equations using traditional finite element techniques 

[5]- 

Traditional methods for determining the radiation heat 
transfer flux q r are based on the assumptions that 
surfaces are isothermal and the incident radiant heat 
fluxes on them are uniform. The isothermal surface 
assumption is inconsistent with finite-element 
formulation since the temperature over the element 
varies according to its shape function. To minimize the 
error introduced by this assumption, the mesh size must 
be controlled to limit the temperature variation over an 
element surface. This can lead to a large number of 
elements for structures with large temperature variation 
throughout. View factors will have to be computed 
between all of these surfaces, a process that requires FT 
computations forN elements. Also, since radiation 
links surfaces throughout a structure, matrices are no 
longer sparse as they are with conduction problems. 
Thus mesh refinement in a radiation heat transfer 
problem can quickly become overwhelming as 
demonstrated in the case study. In addition, Lobo and 
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Emery [6] report that with intense radiation heat fluxes, 
methods employing low-order basis functions can 
produce erroneous results. 

An alternative to mesh refinement (7a -refinement) is to 
increase the order of the basis functions in the finite- 
element formulation (/^-refinement). In a follow-up 
report to their earlier work (cited above), Lobo and 
Emery [7] demonstrate that the errors occurring under 
intense radiation heat flux conditions are due to the 
violation of the Discrete Maximum Principle. They 
show that one way to alleviate this problem is to use 
higher-order basis functions along the surfaces of 
radiating elements. Surana et al [8] have demonstrated 
the effectiveness of /^-version finite element methods 
for conduction problems, but they do not include 
radiation. Kuppurao and Derby [9] use linear and 
quadratic basis functions for pure radiation problems 
arising in crystal growth systems, and develop several 
methods which do not rely upon the isothermal 
assumption. This work takes a slightly different 
approach to implementing higher-order basis functions 
and presents results for test problems representative of 
the wing box analysis. 

RADIATION HEAT TRANSFER WITH 
HIGHER ORDER ELEMENTS 

Fundamental to the accurate implementation of the /?- 
method in heat transfer problems with enclosure 
radiation is the accurate computation of the radiant heat 
flux q r , and if s variation, over a surface. The 
governing equations for the radiant heat flux are of 
integral form. By assuming the surfaces are isothermal 
with uniform radiant heat fluxes, these integral 
equations can be separated into a set of integral 
geometry equations (view factors) and net radiation 
equations (algebraic set of equations). Forgoing these 
assumptions leaves a set of simultaneous integral 
equations that must be solved. Daurelle et al. [10, 1 1] 
shows that solving the set of integral equations leads to 
more accurate results even when using linear basis 
functions for the temperature field. Their work 
indicates that the set of integral equations is slower to 
solve for a given mesh; however, for a given accuracy 
the set of integral equations is twice as fast as the 
traditional approach using the isothermal surface with 
uniform radiant heat flux assumptions. Two methods 
are developed here: the Radiation Sub-Element method, 
which takes advantage of the well-developed traditional 
methods; and the more accurate Integration method, 
which solves a set of integral equations. 


Radiation Sub-Elements (RSE) 

The Radiation Sub-Element approach is the simplest 
and most straightforward technique to implement 
variable surface radiation. Any code that can currently 
handle surface-to-surface radiation can implement this 
procedure without major modifications. This approach 
breaks the radiation elements into sub-elements for the 
surface to surface radiation exchange problem. These 
radiation sub-elements are then treated in the classical 
approach, i.e. assume that they are isothermal and then 
calculate view factors and absorbed heat fluxes. The 
sub-element absorbed heat fluxes are then transferred 
back to the parent element generating a variable 
radiation heat flux along the element. Several different 
methods to approximate the varying heat load on the 
element have been investigated. 

A one-dimensional problem was used to determine how 
accurately the absorbed heat flux should be modeled. 
The energy equation for one-dimensional heat 
conduction in a rod with an applied heat flux cf is: 

kA^ = q'(x)S 
ax 

where 

k is the thermal conductivity, A is the cross sectional 
area of the rod, S is the rod's surface area, and f is the 
variable heat flux expressed as a polynomial in x. 
Integrating twice gives the general solution: 

T(x) = j^\\q p (x)dxdx + c x x + c 2 

The point of interest here is that the temperature 
solution is a polynomial of order p + 2, that is, it is 2 
orders higher than the polynomial describing q . 
Applying this to finite element analysis, given a shape 
function of order m for the temperature distribution 
along a rod element, the radiation heat flux should be 
calculated to an order of m- 2. 

Five methods to transform the discrete heat fluxes on 
the radiation sub-elements to a continuous function on 
their parent element were investigated. Two two- 
dimensional test problems with one-dimensional 
elements were solved using sixth order polynomials for 
temperature interpolation function on the elements. 

The radiation load vectors were calculated to fourth 
order. This required dividing each element into five 
radiation sub-elements yielding five absorbed heat 
fluxes for each parent element. These five heat fluxes 
were then integrated (numerically) with the shape 
functions to generate the radiation load vector fRJ in 
equation 2. 
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The first method models the heat flux along the element 
as a step function. This is consistent with the method 
used for generating the heat fluxes (the classical 
approach assumes the incident radiation to be constant 
along a surface). However, incident radiation does not 
typically vary in this fashion, rather it is usually a 
smooth function which may have significant variation 
depending on the geometry involved. (And if higher 
order elements are being used for the analysis, there 
probably are large variations along the element.) 

The second method approximated the heat flux by 
linearly interpolating between the five heat flux values. 
The heat fluxes were placed at their sub-element 
centroid for the interpolation process. This method 
interpolates between two values when possible, 
otherwise the closest value is used. Values outside of 
two element centroids (i.e., near the open ends of the 
surface) are assigned the same value as the element 
centroid. 

A simple extension of the interpolation method is the 
interpolation/extrapolation method. As the name 
implies, this method simply extrapolates the local heat 
flux value when interpolation is not possible. While 
this method may lead to more accurate solutions in 
many problems, it may also lead to problems if the 
extrapolated values become unrealistic. 

To improve the absorbed radiation heat flux 
approximation a fourth order curve fit was implemented 
next. The five heat fluxes were located at their sub- 
element centroid, and a fourth order polynomial was fit 
to the points. 

The final radiation heat flux approach uses cubic 
splines. Once again the question arises as to how to 
best extrapolate the data when interpolation is not 
possible (i.e. near the end points of an element). Three 
extrapolation methods for the cubic spline were 
investigated. The first approach, also called a natural 
cubic spline, sets the second derivative of the curve to 
zero at the end points. This translates to a linear 
extrapolation of the heat flux values and yielded the 
poorest results. The second method extrapolated the 
data based on a constant second derivative. This 
significantly improved the results. The third method 
allows for a fully cubic curve at the end points and is 
based on a linearly extrapolated second derivative. This 
proved to be the most accurate of all the RSE methods, 
and while no problems were encountered in testing, it 
still extrapolates data and therefore is susceptible to the 
problems associated with such operations. 


Integration Method (IM) 

Regardless of the approach used in the RSE method to 
approximate the radiation heat flux on an element, the 
approximation is based on data from the traditional 
calculation procedures. While this makes 
implementation simple, it still suffers from the 
assumptions used in developing the traditional 
approaches, namely the uniform temperature and 
uniform incident radiation assumptions (albeit they are 
made on sub-elements). 

To study the problem without making the above 
assumptions consider an element k that is part of an 
enclosure of N elements. A heat balance at the point x k 
on element k is: 

q k {x k ) = H k {x k )-J k {x k ) (3) 

where 

H k {x k ) = e k aT k (x k ) + p k J k (x t ) (4) 

dA k J k (x k ) = ± J H j ( Xj , ** )dA, 

M A t 

This last equation can be simplified by applying the 
reciprocity relation 

dAxlF^ = dA k dF dk _ dj 

to give 

l J H / x / F M-d/ x r x k ) 

j 

Substituting this expression into equation 3 gives: 

q k (x k ) = //,(*,)- X J Hj ( Xj )dF^ O j . ) ( 5 ) 

j* 1 Aj 

Equations 3 and 4 can be combined to eliminate the 
surface irradiation term yielding: 

H t (x k ) = oT* (x k ) - — q k (x k ) 

£ k 

which can then be substituted back into equation 5 for 
H k and Hj (with the subscript changed from k to j) to 
give 


° k 

X \ [< (x ) - q (x )lflF (x ,x k ) (6) 

4 £ j 

This equation relates the net surface heat flux to the 
surface temperatures. Note that it is a nonlinear integral 
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equation and that to calculate the heat flux on surface k , 
the surface heat flux on all other surfaces must be 
known in addition to the temperatures of all surfaces 
(including surface k). While this might seem a bit more 
complex, recall that an iterative solution technique is 
generally required to solve the nonlinear finite element 
equations. Thus, the temperatures and heat fluxes from 
the previous iteration values can be used on the right 
hand side of equation 6. 

There are several ways to evaluate equation 6, but 
doing it numerically will only give a single value at 
some point x k . Thus to obtain a function q k (xij, some 
approximation method similar to those shown in the 
RSE method is required. But the radiation heat flux is 
only used to calculate a load vector { R r where it is 
multiplied by the element shape function and integrated 
over the area of the element. If this integration is done 
numerically, the heat flux does not need to be a 
continuous function; it only needs to be evaluated at the 
points defined by the numerical integration scheme. 

The approach taken here is to evaluate the radiation 
load vectors using Gauss Integration. Thus the heat 
flux on an element only need be determined at the 
Gauss points along the element. Also, by using the 
same integration scheme to evaluate the integrals in 
equation 6, no approximating function qk( x k) is 
required. The following steps are used to determine the 
radiation load vector for this solution method: 

1 . Initialize unknowns (temperatures and 
radiation heat fluxes). 

2. Calculate differential form factors 
between Gauss points on all elements. 

(Begin iteration) 

3. Calculate radiation heat fluxes at Gauss 
points on all elements. 

4. Integrate radiation heat fluxes to obtain 
radiation heating load vector. 

5. Solve for and update temperatures. 

6. Check convergence and repeat iteration if 
necessary (go to step 3). 

APPLICATION 

Two-Dimensional Test Problems 

To investigate these techniques, two simple test 
problems were developed. The test problems involve 
two one-dimensional elements in radiative equilibrium. 
The first test problem, shown in Figure 6, represents 
heat transfer in the comer of a rectangular enclosure. 


The two elements are oriented at right angles with each 
other and share a common end point location (there is a 
separate node at this location for each element). The 
vertical element is held at a constant temperature of 
1000°R and the horizontal element is allowed to come 
to radiative equilibrium. The horizontal element’s 
length was set to 5 times the vertical element length in 
order to allow a reasonable temperature gradient to 
develop along the element. Both surfaces were 
modeled as black bodies so that an exact solution could 
be found. 



1 000°R 


5L 


x 


Figure 6: Schematic of test problem 1. 


The second test problem is shown in Figure 7. This 
problem considers two parallel black plates each of 
length L separated a distance d apart. Once again the 
top plate was held at a constant temperature and the 
other plate allowed to come to radiative equilibrium. 
Results presented here are for a d/L ratio of 0.1. 


T = 1 000°R 



d 


L 

Figure 7: Schematic of test problem 2. 

Results 

Radiation Sub-element Method 

Results for test problem 1 using the RSE methods are 
shown in Figure 8. Results for all the methods of 
approximating the radiation heat loads are shown along 
with the exact solution. Although the step function 
method is consistent with the isothermal/uniform heat 
flux assumption used to generate the heat flux data, the 
resulting temperature prediction was poor with the error 
approaching 20% at the element end points. 

Results of the interpolation method are slightly better 
than the step function results; however, both methods 
significantly under-predict the peak temperature at x/L 
= 0. The problem here is that the incident heat flux is 
rising rapidly (x/L = 0 corresponds to the comer 
location). Both the step function and linear 
interpolation methods model the heat flux as constant in 
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this area (step function holds the value constant from 
x/L = 0 to 0.2, linear interpolation function holds the 
value constant from x/L = 0 to 0.1) and thus under- 
predict the actual temperature. Extrapolating the data 
near the endpoints improves the results considerably 
cutting the error at these locations by approximately 
50%. The extrapolated heat flux values still under 
predict the actual heat flux values near x/L=0, and over 
predict the actual heat flux values near x/L=5. 



Figure 8: Results for test problem 1. 


The fourth order curve fit gives excellent results over 
most of the element; however, there is a significant 
error in the temperature at x/L of 5. This behavior is 
due to the nature of the fourth order polynomial defined 
by the 5 heat flux values. While the actual radiation 
heat flux for this problem continues to decrease slowly 
as x/L approaches 5, the polynomial fit begins to 
increase here. Several approaches could be 
implemented to avoid this problem; however, while 
extrapolating the results is causing errors at x/L=5, it is 
helping yield a better answer at x/L=0. One approach 
that would certainly improve the results for this case 
would be to subdivide the element into additional 
radiation sub-elements and use a curve fit to smooth out 
the polynomial. The problem here is the additional 
computations required for the new sub-elements. 

Finally, the cubic spline curve fit gave the best overall 
results. The results shown in Figure 8 are for the fully 
cubic spline. The temperature results show good 
agreement over the whole element with the largest 
errors occurring near x/L=5. 

In general all the approximation methods suffer near the 
endpoints, regions where the data must be extrapolated. 
These methods may be improved by calculating the 
heat flux at the endpoints of the element. This would 
require computation of radiation heat fluxes at points, 
as opposed to finite surfaces. The traditional methods 
used to calculate the heat fluxes on an element would 
have to be modified to handle both finite areas as well 
as points; however, this complication is only minor and 
may well be worth the additional effort. Another 


consequence which might prove more significant is that 
the interior points would be farther apart (assuming 
only 5 points) possibly increasing the approximation 
error in the interior of the element. 



Figure 9: Results for test problem 2. 


Figure 9 shows the results of these methods applied to 
test problem 2 along with the exact solution. The 
results are similar to those from test problem 1 with the 
cubic spline curve fit method giving the best overall 
results. To quantify the overall performance of the 
various methods, the following error indicator was 
calculated for each method: 


error indicator = 


T.{Tn-Tjf 


V XX, 


where T FE is the temperature from finite element 
calculation, is the temperature from the analytical 
solution, and i indicates the number of the & data point 
shown in the figures (total of 21 points) 

This error indicator is shown in Figure 10 for test 
problem 1 and in Figure 1 1 for test problem 2. As 
expected the cubic spline curve fit has the lowest error 
indicator value. 



step function linear interp. interp/extrap 4th order of cubic spline 


Figure 10: Error indicator for test problem 1. 
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step function linear interp. interp/extrap 4th order cf cubic spline 


Figure ll:Error indicator for test problem 2. 

Integration Method 

The results for test problem 1 using the integration 
method are shown in Figure 12. The figure shows 
results using a 7 Gauss point integration scheme, which 
performs so well that it is difficult to distinguish 
between the integration method and the exact solution. 
The results for the cubic spline radiation sub-element 
method are also included to gage how well the 
integration method performs against the best radiation 
sub-element method. The integration method has 
nearly an order of magnitude reduction in the error 
indicator. 



Results for test problem 2 are shown in Figure 13. 

Once again results for the cubic spine radiation sub- 
element method are included for comparison. Note that 
the integration method slightly over-predicts the 
temperature at the center of the element while the 
radiation sub-element method slightly under-predicts 
the temperature there. The integration method used 16 
Gauss points to obtain this solution, significantly more 
than was necessary for the first test problem. This 
higher accuracy scheme was necessary because of the 
geometry associated with test problem 2, namely large 
surfaces separated by a small distance. This is a 


common problem in radiation heat transfer and view 
factor calculations and is discussed in many view factor 
references. 



Figure 13: IM vs. RSE results fro test problem 2. 

Several additional finite element analyses were 
performed to compare these methods to traditional 
methods for computing radiation heat transfer. Each 
test problem was analyzed using the traditional finite 
element approach-linear shape functions for the 
temperature field with the radiation heat flux assumed 
constant over each element. Each surface was 
subdivided into a number of linear elements and the 
mesh was refined until the error indicator was reduced 
to the level produced by the higher-order methods. 

Table 1 lists the number of linear elements required to 
match the solution accuracy of the new methods. The 
traditional approach required 10 elements for test 
problem 1 and 1 7 elements for test problem 2 to match 
the error indicator for a single element using the cubic 
spline radiation sub-element method. Similarly, the 
traditional approach required 46 elements for test 
problem 1 and 42 elements for test problem 2 to match 
the error indicator for a single element using the 
integration method. These traditional analyses were 
carried out using a totally separate computer code 
previously developed by the author. This code had very 
little similarity to the code used to develop the higher 
order methods, so no attempt was made to compare the 
computational costs (run times) for the methods. 

Table 1: Number of linear elements required to 
match the accuracy of a single p = 6 element using 
the RSE and IM methods. 


Test 

problem 

Radiation sub-element 
method 

(cubic spline curve fit) 

Integration 

method 

1 

10 

46 

2 

17 

42 
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Three-Dimensional Test Problem 

Based on the results from the two-dimensional test 
problems, the cubic spline curve fit RSE and integration 
methods were implemented in a three-dimensional 
finite element code. The test problem for this case, 
shown in Figure 14, consists of two parallel six-inch 
square plates separated by a distance of six inches. The 
top plate is held at 1000°F and there is no conduction in 
either plate. The emissivity of both plates is 1.0 so that 
the analytical solution could be obtained: 


r(x,y) = 1000 


36 


;r[(x-x) 2 + {y-y) 2 +36} 


-dx dy 


A surface plot of this solution is shown in Figure 15. 

A slightly different method was used to measure the 
accuracy of the solutions for this test problem so that 
comparison with ^-refinement could be made. The L 2 
norm of the error for element K is: 


ill'll L 2 {K) 1/C analytic '^'finite element ) dx dy 

The error for the problem is found by taking the square 
root of the sum of all the element errors in the domain 
Q: 

Solutions for this problem were generated using h- 
refmement with the traditional radiation approach 
(isothermal surfaces with uniform radiation heat flux); 
L 2 norms for these cases are shown in Table 2. The 
amount of mesh refinement was limited by limitations 
in the view factor software; however, the cases 
presented give a good indication of how quickly the 
error drops with /^-refinement. Results for the radiation 
sub-element method are presented in Table 3. The data 
indicate that increasing the interpolation function gives 
more accurate solutions for a given number of degrees 
of freedom. However, the integration method provides 
the best accuracy as the data in Table 4 show. Only the 
even polynomial results are presented because the 
solution is an even function, and the odd polynomials 
do not improve the solution accuracy. Note that using a 
single element with a polynomial of order 2 (in each 
direction) with the integration method provides a 
solution more accurate than any of the other 
approaches. For nine degrees of freedom, the h~ 
refinement L 2 error norm is 19.45, the radiation sub- 
element method with nine degrees of freedom has an L 2 
error norm of 1.97 and the integration method with nine 
degrees of freedom has an L 2 error norm of 0.23. In all 
cases, for a given number of degrees of freedom, the 
integration method produces the lowest error followed 


by the radiation sub-element method with the h- 
refinement method giving the largest error. 



Figure 14: Parallel plates under radiation 
equilibrium conditions with the lower plate held at a 
uniform temperature. 



Figure 15: Analytical solution for 3-D test problem. 

Table 2: Error norms for h-refinement using the 
traditional radiation method. 


Number of 
p=l elements 

Number of 
degrees of 
freedom 

L 2 error 
norm 

1 

4 

19.4522 

4 

9 

19.4522 

9 

16 

9.74385 

16 

25 

6.51942 

25 

36 

4.65848 

36 

49 

3.56349 

49 

64 

2.83546 

64 

81 

2.32688 

81 

100 

1.95422 

100 

121 

1.67149 

400 

441 

0.59645 
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Table 3: Error norms for /7-refinement using the 
radiation sub-element method. 


Order of 
element ( p ) 

Number of 
degrees of 
freedom 

L 2 error 
norm 

1 

4 

19.4522 

2 

9 

1.96614 

3 

16 

1.10128 

4 

25 

0.73404 

5 

36 

0.52810 

6 

49 

0.40612 


Table 4: Error norms for /7-refinement using the 
integration method. 


Order of 
element {p ) 

Number of 
degrees 
Of freedom 

L 2 error 
norm 

1 

4 

19.4547 

2 

9 

0.23185 

4 

25 

0.00789 

6 

49 

0.00081 


CONCLUSIONS 

Two methods for determining internal radiation heat 
transfer have been developed for higher-order finite 
elements. The first method divides the higher-order 
element into a number of sub-elements, calculates the 
radiant heat flux on the sub-elements using traditional 
methods, and then curve fits this data to determine the 
radiation heat flux along the higher-order parent 
element. The second method numerically integrates the 
radiation heat transfer equations over the higher-order 
element using an efficient Gaussian integration scheme. 
Comparisons with analytical solutions show that the 
integration scheme is generally more accurate than the 
sub-element method. Comparisons of these results to 
those of traditional linear finite elements demonstrate 
the potential for improved computational performance 
given a required level of accuracy. 
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